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ABSTRACT 

Aims. Using a reconstruction of sunspot numbers stretching over multiple millennia, we analyze the statistics of the 
occurrence of grand minima and maxima and set new observational constraints on long-term solar and stellar dynamo 
models. 

Methods. We present an updated reconstruction of sunspot number over multiple millennia, from 14 C data by means 
of a physics-based model, using an updated model of the evolution of the solar open magnetic flux. A list of grand 
minima and maxima of solar activity is presented for the Holocene (since 9500 BC) and the statistics of both the length 
of individual events as well as the waiting time between them are analyzed. 

Results. The occurrence of grand minima/maxima is driven not by long-term cyclic variability, but by a stochas- 
tic/chaotic process. The waiting time distribution of the occurrence of grand minima/maxima deviates from an ex- 
ponential distribution, implying that these events tend to cluster together with long event-free periods between the 
clusters. Two different types of grand minima are observed: short (30-90 years) minima of Maunder type and long 
(>110 years) minima of Sporer type, implying that a deterministic behaviour of the dynamo during a grand minimum 
defines its length. The duration of grand maxima follows an exponential distribution, suggesting that the duration of 
a grand maximum is determined by a random process. 

Conclusions. These results set new observational constraints upon the long-term behaviour of the solar dynamo. 
Key words, long-term solar activity - cosmic rays - solar dynamo - sunspot numbers 



1. Introduction 

The Sun is the only star whose magnetic activity can 
, be studied on long time scales. Direct solar observations 
since 1610 reveal great variability of the cycle averaged 
magnetic activity level of the Sun - from the extremely 
quiet Maunder minimum (second half of 17th century) 
" up to the modern episode of enhanced activity since 
' the middle of the 20th century. The Maunder mini- 
mum is representative of grand minima of solar activity 
(e.g., |Eddy 1977afl , when sunspots almost completely 
vanished from the solar surface, while the solar wind 
appear ed to continue b l owing, although at a reduced 
pace (|Cliver et al. 19981 lUsoskin et al. 200ip . A grand 
minimum is believ ed to corresp o nd to a special stat e 
of the dynamo (jSokoloff 20041 |Miyahara et al. 2006[ ), 
and its very existence poses a challenge for solar dy- 
namo theory. It is noteworthy that dynamo models do 
not agree how often such episodes occur in the Sun's 
history and whether their appearance is regular or 
random. For example, the commonly used mean-field 
dynamo yields a fairly regular 11-year cycle, while 
there are also d ynamo models in cluding a stochas- 
tic driver (e.g., IChoudhuri 19921 ISchlnitt et al. 19961 

jWeiss fc Tobias 20001 

ICharbonneau 20011 2004) which predict 
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intermittency of the solar magnetic activity. The presence 
of grand maxima of solar activity has been mentioned 
QEddy 1977a] lUsoskin et al. 20031 ISolanki et al. 2004|) but 
has not yet been studied in great detail. 

Thanks to the recent development of precise technolo- 
gies, including accelerator mass spectrometry, solar activity 
can be reconstructed over multiple millenia from concen- 
trations of cosmogenic isotopes 14 C and 10 Be in terrestrial 
archives. This allows one to study the temporal evolution 
of solar magnetic activity, and thus of the solar dynamo, on 
much longer time scales than available from direct measure- 
ments. Consequently, a number of attempts to investigate 
the occurrence of grand minima in the past, using radio- 
carbon 14 C data in tree rings, have been undertaken. E.g., 
Eddy (1977a) identified major excursions in the available 
14 C record as grand minima and maxima of solar activity 
and presented a list of 6 grand minima and 5 grand max- 
ima for the last 5000 years. Stuiver & Braziunas (1989) 
studied grand minima as systematic excesses of the high- 
pass filtered 14 C record and suggested two distinct types 
of the grand minima: shorter Maunder-type and longer 
Sporer- like minima (cf. IStuiver et al. 1991]) . Later Voss et 
al. (1996) defined grand minima in a similar manner and 
provided a list of 29 such events for the last 8000 years. A 
similar analysis of excursions of the 14 C production rate has 
been presented by Goslar (2003). However, because of the 
lack of adequate physical models relating the radicarbon 
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abundance to the solar activity level, such studies retained 
a qualitative element. The use of high-pass filtered 14 C data 
is based on the assumption that solar activity variations 
are important only on short times, while all the long-term 
changes in radiocarbon production are attributed solely to 
the slowly changing geomagnetic field. This method ignores 
any possible long-term changes in the solar activity (e.g., 
on time scales longer than 500 years for Voss et al. 1996). 
There is, however, increasing evidence that solar activity 
varies on multi-centennial to multi-millennial time scales 
(jMcCracken et al. 20041 lUsoskin et al. 2006a|) . A recently 
developed approach, based on physics-based modelling of 
all links relating the measured cosmogenic isotope abun- 
dance to the level of solar activity, allows for quantitative 
reconstruction of the solar activity level in the past, and 
thus, for a more realistic definition of the periods of grand 
minima or maxima. 

Here we study the statistics of occurrence of grand min- 
ima/maxima throughout the Holocene and impose addi- 
tional observational constraints on the dynamo models of 
the Sun and Sun-like stars. 

2. Past solar activity 

Solar activity on multi-millcnial time scales has been re- 
cently reconstructed using a physics-based model from 
measurements of 14 C in tree rings (see full details in 
ISolanki et al. 2004[ lUsoskin et al. 2006ap . The validity of 
the model results for the last centennia has been proven 
by independent data on measurements of 44 Ti in stony 
meteorites ()Usoskin et al. 2006b|) . The reconstruction de- 
pends on the knowledge of temporal changes of the geo- 
magnetic dipole field, which must be estimated indepen- 
dently by paleomagnetic methods. Here we compare two 
solar activity reconstructions, which are based on alterna- 
tive paleomagnetic models: one which yields an estimate of 
the virtual aligned dipole moment (VADM) since 9500 BC 
( |Yang et al. 20 00), and the other a recent paleomagnetic 
reconstruction of the true dipole moment since 5000 BC 
(|Korte fc Cons table 2 005|) . We note that the geomagnetic 
dipole moment obtained by Korte & Constable (2005) lies 
systematically lower than that of Yang et al. (2000), lead- 
ing to a systematically higher solar activity reconstruction 
in the past (jUsoskin et al. 2006a[) . While the geomagnetic 
reconstruction of the VADM by Yang et al. (2000) provides 
an upper bound for the true dipole moment, the more re- 
cent work of Korte & Constable (2005) may underestimate 
it. Thus we consider both models as they bound a real- 
istic case. We note that the Yang et al. (2000) data run 
more than 4000 years longer and give a more conservative 
estimate of the grand maxima. 

Different indices of the reconstructed solar activity 
are shown in Fig. [1] Most directly related to the 14 C 
production in the atmosphere is the modulation poten- 
tial cj) (e.g., |Castagnoli fe Lai 1980[ IMasarik fc Beer 19991 
lUsoskin et al. 2002|) whose variations are shown in panel 

A. The modulation potential is a parameter describing 
the spectrum of galactic cosmic rays (see definition and full 
description of this index in lUsoskin et a l. 2005 ) . Using a 
model of the heliospheric transport of cosmic rays, the mod- 
ulation potential can be nearly linearly related to the open 
solar magnetic flux F a . The reconstructed long-term vari- 
ations of the open solar magnetic flux are shown in panel 

B. Next, using a model of the open magnetic flux forma- 
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Fig. 1. Long-term solar activity reconstruction from 
14 C data. All data are decadal averages. Solid (de- 
noted as 'Y00') and grey ('K05') curves are based 
on the paleo-geomagnetic reconstructions of Yang et 
al. (2000) and Korte & Constable (2005), respec- 
tively. A) The modulation potential <f>. Big circles 
('NM') denote the <j) values obtained from direct cos- 
mic ray measurements since 1951 (lUsoskin et al. 2 005). 

B) Open flux F . Reconstruction from sunspot numbers 
QKrivova, Balmaceda fc Solanki 2007| ) is used after 1610. 

C) Sunspot numbers reconstructed here. The Y00 and 
K05 curves are called SN-L and SN-S series, respectively 
throughout the paper. Observed group sunspot numbers 
dHoyt fc Schatten 1998[ ) are shown after 1610. 



tion, one can estimate the sunspot numbers from the F a 
data. This was done earlier using the open flux model by 
Solanki et al. (2000, 2002). However, the model relating 
sunspots to the open magnetic flux, has been updated re- 
cently by Krivova, Balmaceda & Solanki (2007). Starting 
from sunspot numbers (SN) they computed open and to- 
tal magnetic flux as well as the total solar irradiance (from 
SN and total magnetic flux). Besides the active regions, the 
influence of ephemeral active regions is included. They de- 
termined the free parameters of the models used by requir- 
ing the model output to reproduce the best available data 
sets (of open and total fluxes as well as total solar irradi- 
ance) with the help of a genetic algorithm. In particular, 
the improved total flux data set of Arge et al. (2002) and 
the total irradiance composite of Frohlich (2006) set tight 
constraints, which could only be met by revising the rela- 
tionship between SN and total emergent magnetic flux. This 
revised relationship is employed here too. Accordingly, here 
we use this updated model to convert F Q into the sunspot 
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Fig. 2. Scatter plot of decadal sunspot numbers for 9500 
BC - 1900 AD published by Solanki et al. (2004), i? G (S04), 
vs. sunspot numbers obtained in this work in a way identi- 
cal to S04 but using the updated open solar flux model of 
Krivova, Balmaceda & Solanki (2007). 



number (Eq. IA.17)) . which is described in detail in Appendix 
A. 

A comparison between the sunspot number series ob- 
tained using the new open flux model with those published 
by Solanki et al. (2004) is shown in Fig. [2] The scatter in 
the relationship is very small (correlation coefficient 0.995), 
but there is a small systematic difference between the two 
series. The newly obtained sunspot numbers are slightly 
higher, with the mean difference being about 1.6 and the 
standard deviation about 2. At large SN values, however, 
the opposite is observed, with the new reconstructed SN 
being smaller by values up to 14. This reduces the largest 
peaks (grand maxima) somewhat. This difference is a result 
of the different relationship between SN and total emerging 
magnetic flux in active regions used by Krivova, Balmaceda 
& Solanki (2007) than by Solanki et al. (2000, 2002). We 
note that the difference is kept within the uncertainty of 
the reconstruction, which is about 8 for the last millenia 
and up to 15-20 in the be ginning of the Hol ocene (see 
Supplementary Material to ISolanki et al. 2004]) . The dif- 
ference between the original sunspot numbers published in 
Usoskin et al. (2006a) based on the paleomagnetic model 
by Korte & Constable (2005) and those obtained here us- 
ing the updated open flux model has a mean of 1.0 and 
standard deviation about 2.4. 

The 11,000-yr long data sets of the decadal sunspot 
number similar to that by Solanki et al. (2004) but with the 
updated open flux model is shown in Fig. \Vp. It is called 
the SN-L series throughout the paper. The shorter series 
(Fig. |Tp) , which is similar to that by Usoskin et al. (2006), 
is called SN-S henceforth. Afte r 1610 AD, the actually ob- 
served group sunspot numbers ( Hoyt fc Schatten 1998[ ) has 
been used instead of the reconstructions. 




-4000 -3000 -2000 -1000 1000 2000 

Years (-BC/AD) 

Fig. 3. Sunspot activity SN-L throughout the Holocene 
(see text) smoothed with a 1-2-2-2-1 filter. Blue and red 
areas denote grand minima and maxima, respectively. The 
entire series is spread over two panels for better visibility. 

Before identifying the grand minima and max- 
ima, the decadal resolution data have been smoothed 
with the Gleissberg (1-2-2-2-1) filter, which is regu- 
larly applied when studying long-term variations of 
solar activity in order to suppress the noise (e.g. 



Gleissberg 1944 


|Soon, Posmentier & Baliunas 1996} 


Mursula & Ulic 


i 


1998). In order to check the effect of 



the filter we have studied a number of artificial SN series 
containing a total of 1000 grand minima of 60-yr duration 
(at the level of SN<15) each. A noise with a = 10 has 
been added to the series, and the grand minima have been 
identified again as SN<15 in both the raw and 1-2-2-2-1 
smoothed noised series. We found that 35% of grand 
minima are incorrectly identified (too short, too long 
or split in two short episodes, comparing to the "real" 
signal) in the raw noisy series. The filtering reduces the 
mis-identification to 13±3%, i.e. 3-fold. Thus, the use of 
the 1-2-2-2-1 Gleissberg filter reduces the effect of noise 
on the grand minima/maxima definition and makes the 
results more robust. This smoothing, however, leads to 
a reduction in the amplitude and a slight underestimate 
(about 7% according to the above numerical experiment) 
of the number and duration of short, less than 30 years 
long minima and maxima. 

The filtered SN-L and SN-S series are shown in Figs. [3] 
anddJ respectively. We analyze both reconstructed SN data 
sets in equal details. 

3. Definitions 

3.1. Grand minima 

We have defined a grand minimum as a period when the 
(smoothed) SN level is less than 15 during at least two 
consecutive decades - this corresponds to blue-filled areas 
in Figs. [3] and [4] However, taking into account the uncer- 
tainty of the SN reconstruction and the influence of the 
filtering, we have also considered as minima clear dips in 
the SN whose bottom is between 15 and 20 and the depth 
(with respect to surrounding plateaus/maxima) exceeds 20. 
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Fig. 4. Sunspot activity SN-S (see text) smoothed with a 
1-2-2-2-1 filter. Blue and red areas denote grand minima 
and maxima, respectively. 

Therefore, such periods as, e.g., ca. 6400 BC are considered 
as grand minima (see Fig. [3]) even though their bottoms 
are above 15. On the other hand, the period ca. 4450 BC 
is not counted as a grand minimum because its minimum 
lies above 15 and its depth is less than 20. All 27 grand 
minimum periods thus defined in the SN-L series are listed 
in Table [1] together with their approximate duration, de- 
fined as the period of time when SN was below 15 (20) 
as discussed above. Among them there are two periods (ca 
9200 BC and 7500 BC) when the dominant grand minima 
are interrupted by a 1-2 decades long upward excursions. 
We regard these period as continuous Sporer-like minima. 
Together these grand minima have a total duration of 1880 
years, so that the Sun spends about 17% of the time in a 
grand minimum state. 

We note that all the grand minima after 3000 BC dis- 
cussed by Eddy (1977a, 1977b) are present in TableHJ which 
however contains also minima at ca. 1040 BC and 2860 BC 
not found by Eddy. On the other hand, all the grand min- 
ima listed in the Table are mentioned by Voss et al. (1996), 
but the latteiQ list more minima, e.g., three minima be- 
tween 200 BC and 200 AD which do not appear in our 
series. Most of the minima listed here can also bee seen in 
Fig. 2 of Goslar (2003). Therefore, we conclude that our 
definition of grand minima applied to the present data set 
gives results generally in agreement with earlier studies but 
not identical to them. In particular, it gives more details 
than the study by Eddy (1977a) but discards some small 
fluctuations mentioned by Voss et al. (1996). 

The grand minima listed in Table Q] dating after 5000 
BC are identical for both SN-L and SN-S series (except for 
one minimum ca. 1385 BC in the SN-S series which does 
not match the formal definition). These grand minima will 
be used for the further analysis. 

3.2. Grand maxima 

Similar to Solanki et al. (2004), we define as a grand max- 
imum of solar activity a period when SN exceeds 50 dur- 
ing at least two consecutive decades (see red filled areas in 
Figs. [3] and[4|). If two consecutive maxima are separated 
by less than 30 years they are considered as a single max- 
imum (e.g., ca. 9000 BC in the SN-L series), i.e. they are 
treated in a way similar to grand minima. We have iden- 
tified 19 grand maxima (of a total duration of 1030 years, 
corresponding to about 9% of the time) in the SN-L series 

1 Definition of minima by Voss et al. (1996) was based solely 
on the relative variations of A 14 C. 



Table 1. Approximate dates (in -BC/AD) of grand minima 
in the SN-L series (see text). 



No. center duration comment 



1 


1680 


80 


Maunder 


2 


1470 


160 


Sporer 


3 


1305 


70 


Wolf 


4 


1040 


60 


a) 


5 


685 


70 


b) 


6 


-360 


60 


a,b,c) 


7 


-765 


90 


a,b,c) 


8 


-1390 


40 


b,e) 


9 


-2860 


60 


a,c) 


10 


-3335 


70 


a,b,c) 


11 


-3500 


40 


a,b,c) 


12 


-3625 


50 


a,b) 


13 


-3940 


60 


a,c) 


14 


-4225 


30 


c) 


15 


-4325 


50 


a,c) 


16 


-5260 


140 


a,b) 


17 


-5460 


60 


c) 


18 


-5620 


40 




19 


-5710 


20 


c) 


20 


-5985 


30 


a,c) 


21 


-6215 


30 


c,d) 


22 


-6400 


80 


a,c,d) 


23 


-7035 


50 


a,c) 


24 


-7305 


30 


c) 


25 


-7515 


150 


a,c) 


26 


-8215 


110 




27 


-9165 


150 





a) Discussed in Stuiver (1980) and Stuiver & Braziunas (1989). 

b) Discussed in Eddy (1977a, 1977b). 

c) Shown in Goslar (2003). 

d) Exact duration is uncertain. 

e) Does not appear in the SN-S series. 

since 9500 BC, including also the modern maximum. These 
are listed in Table [2j Four out of six grand maxima found 
here after 3000 BC coincide with those pointed out by Eddy 
(1977a, 1977b). In the SN-S series, 23 grand maxima (of a 
total duration of 1560 years, corresponding to about 22% of 
the time) are identified since 5000 BC, as listed in Table [3l 
All maxima identified in the SN-L series are present also in 
the SN-S series, but the latter yields more maxima satisfy- 
ing the same definition before 1500 BC (after ca. 1500 BC 
the maxima are nearly identical). This indicates that the 
identification of maxima is less robust than for grand min- 
ima, and is more dependent on the definitions and model 
assumptions. 

3.3. Waiting time distribution 

The interval between two consequent events is called the 
waiting time. The statistical distribution of waiting times 
(WTD - waiting time distribution) reflects the nature of 
a process which produces the studied events. For instance, 
an exponential WTD is a clear indicator of a purely ran- 
dom, "memoryless" process (e.g., Poisson process), when 
the behaviour of a system does not depend on its preceding 
evolution on both short or long time-scales. Any significant 
deviation of the WTD from an exponential law implies that 
the probability of an event to occur is not time-independent 
but is related to the previous history of the system. We note 
that the occurrence of events generally is random also for 
a non-exponential distribution, but the probability is not 
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Table 2. Approximate dates (in -BC/AD) of grand max- 
ima in the SN-L series. 



No. 


center 


duration 


comment 


1 


1960 


80 


modern, b) 


2 


-445 


40 


- 


3 


-1790 


20 


a) 


4 


-2070 


40 


- 


5 


-2240 


20 


a) 


6 


-2520 


20 


a) 


7 


-3145 


30 


- 


8 


-6125 


20 


- 


9 


-6530 


20 


- 


10 


-6740 


100 


- 


11 


-6865 


50 




12 


-7215 


30 




13 


-7660 


80 




14 


-7780 


20 




15 


-7850 


20 




16 


-8030 


50 




17 


-8350 


70 




18 


-8915 


190 




19 


-9375 


130 




a) Discussed in 


Eddy (1977a, 1977b)). 


b) Cent 


er and duration of the modern maxi 


since it is still ongoing. 




Table 3. Approximate dates (in -BC/ 


ima in 


the SN-S series. 




No. 


center 


duration 


comment 


1 


1960 


60 


a,b) 


2 


-265 


70 




3 


-455 


70 


a) 


4 


-595 


90 




5 


-1065 


50 




6 


-1560 


60 




7 


-1640 


40 




8 


-1775 


170 


a) 


9 


-1995 


210 


a) 


10 


-2165 


30 




11 


-2235 


50 


a) 


12 


-2350 


80 




13 


-2515 


50 


a) 


14 


-2645 


30 




15 


-2715 


50 




16 


-2790 


40 




17 


-2960 


60 




18 


-3030 


40 




19 


-3170 


120 


a) 


20 


-3415 


50 


a) 


21 


-3840 


60 




22 


-4090 


60 




23 


-4870 


20 





a) Exists also in the SN-L series (Table [2}. 

b) Center and duration of the modern maximum are preliminary 
since it is still ongoing. 



uniform in time. This can be interpreted in different ways: 
self-organized criticality (e.g., Ide Carvalho fc Prado 2000, 
IFreem an et al . 20001) . time-dependent Poisson process (e.g., 
[Wheatland 2 003]), some memory in the driving process 
(e.g., |Lepreti et al. 200H|Mega et al. 20030 . The most typ- 
ical non-exponential WTD is a power law which is, e.g., 
a necessary but not sufficient indication of self-organized 
criticality (|de Carvalho fc Prado 20000 . A power law im- 



plies higher tails of the distribution, i.e. higher probabil- 
ity (relative to the exponential function) of occurrence of 
both long and short intervals between the events. A power 
law distribution of the waiting time has been obtained 
for many solar and terrestrial indices on different time 
scales from minutes to 10 5 years: e.g., intervals between 
major earthquakes (jBak et al. 2 002. Me ga et al. 2003 ); in- 
tervals between suc c essive solar flares (jPearce et al~ 1993, 
IBoffetta et al. 19991 IMoon et al. 20010 : waiting time be- 
tween successive coro nal mass ejections ([Wheatland 2003, 
Bcrh ondo et a l. 2006); intervals between bursts in the so- 
lar wind (Freema n et al. 2000 ): repetition time of geomag- 
netic disturbances ( |Papa et al. 20060; intervals betwee n the 
geomagnetic field reversals ( Ponte-Neto fc Papa 20060 , etc. 
Note that many of these processes, which depict different 
degrees of self-organization, are related to energy accumu- 
lation and release. 

Here we study the WTD of the occurrence of grand 
minima and maxima of solar activity in order to understand 
the nature of its long-term evolution. The waiting time is 
defined as the length x of an interval between centers of 
consequent events. 

First we studied the differential distribution which is 
defined as 



N{xi,x 2 } 



(1) 



where N{x\ 1 X2\ is the number of events with the wait- 
ing time x± < x < X2- Statistical errors of the differential 
distribution are estimated from the Poisson statistics as 

VN/(X 2 -X 1 ). 

Since the statistics are poor (19-27 events) and differ- 
ential WTD histograms are rough, we have studied also the 
normalized cumulative distribution defined as 



Y(x) 



N{x,oo} 
7V{0,oo}' 



(2) 



which corresponds to the probability of finding a waiting 
time exceeding x. In this case the statistical errors cannot 
be determined since the points of the cumulative distribu- 
tion are not independent. 

The exponential WTD model is defined as 



y{x) cx exp 



i ; Y(x) cx exp ^ 



T. 



The power law WTD model is defined as 
y(x) cx a; -7 ; Y(x) cx x~ r ; 7 = T + 1. 



(3) 



(4) 



The power law serves mainly to emphasize deviations from 
a purely exponential distribution, since the poor statistics 
hardly allow us to distinguish between a power law and 
other more-stretched-than-exponential distributions, e.g., 
log-Poisson. 

We note that short intervals (shorter than a century) 
cannot be reliably defined because of noise and filtering. 
Statistics of very long intervals is not reliable either because 
of the limited length of the analyzed series. Therefore, when 
fitting the data we will ignore the shortest and longest inter- 
vals, i.e. first and last points of the cumulative distribution 
(the number of bins of the differential distribution is left 
unchanged) . 
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Fig. 5. Histogram of the sunspot number SN series re- 
constructed here for 9,500 BC - 2000 AD. Hatched areas 
correspond to directly observed sunspots after 1610. The 
curve represents the best fit normal distribution. 

4. Analysis and results 

4.1. Sunspot number distribution 

First we have constructed histograms of the sunspot num- 
bers. The histogram for the SN-L series is shown in Fig. [5] 
While being close to a normal distribution (mean=31, 
a = 30), there is an apparent excess both at very low 
sunspot numbers, corresponding to the grand minima, and 
at very high values, corresponding to grand maxima. The 
overall distribution is consistent with the direct observa- 
tional record after 1610, suggesting that the latter is a rep- 
resentative sample for the sunspot activity statistics, in- 
cluding grand minimum and maximum 0. This distribution 
with these excesses suggests that grand minima and max- 
ima are special states of the solar dynamo that cannot be 
explained by random fluctuations or noise in the data (see 
also forthcoming sections). 

4.2. Grand minima 

4.2.1. Waiting time distribution 

The distribution of the waiting time between grand minima 
is shown in the left-hand panel of Fig. [6] together with the 
best fit approximations. Parameters of the best-fit approxi- 
mations are shown in TablelU (row A). The best fit exponen- 
tial model (Eq. [3]) gives t = 330 ± 50 years, which roughly 
corresponds to the mean frequency of grand minima occur- 
rence. The exponential model agrees only relatively poorly 
with the observed WTD. The best fit power law model 
(Eq. 2]) agrees reasonably with the observed WTD. 

The cumulative WTD is shown in the right-hand panel 
of Fig. [6]together with the best fit approximations (TablelU 
row A, columns 4 and 5). The power law model agrees 
well with the bulk of the data except for the very far tail 
(x >1000 years). However, this tail contains only two events 
and is not representative. (The x 2_ test cannot be applied 
to the cumulative distribution since the points are not in- 
dependent.) The exponential model poorly reproduces the 

2 Note that in Fig. \5\ intermediate SN values are seemingly 
underepresented in modern times, but this is due to the loga- 
rithmic scale. 




Interval length (Years) Interval length (Years) 

Fig. 6. Differential (left panel) and cumulative (right 
panel) distribution of the waiting time between subsequent 
grand minima. The histogram (left) and circles (right) rep- 
resent the observed distribution, while solid and dotted 
lines depict best fit power law and exponential approxi- 
mations, respectively. 



WTD. As an additional test we compare the parameters 
of the models describing differential and cumulative distri- 
butions, viz. T« t and 7 ~ T + 1. Both models pass this 
test. 

We conclude that a power law model better describes 
the observed WTD for grand minima, although an expo- 
nential decay cannot be completely ruled out. This is valid 
also for the SN-S series whose grand minima (except of one 
ca. 1385 BC) coincide with those in the SN-L series after 
5000 BC. 

4.2.2. Duration of grand minima 

A histogram of the duration of grand minima (Table [I) is 
shown in Fig. [7J The mean duration is 70 year but the 
distribution is not uniform. The minima tend to be ei- 
ther of a short duration, between 30 and 90 years simi- 
lar to the Maunder minimum, or rather long, longer than 
110 years similar to the Sporer minimum. This agrees with 
the earlier conclusion on two different types of grand min- 
ima (jStuiver fc Braziunas 19891 IGoslar 2003) . This sug- 
gests that a grand minimum is a special state of the dynamo 
whose duration is not random but is defined by some intrin- 
sic process. Note, however, that only 3 of the 5 Sporer-likc 
minima are clear long grand minima while the other 2 are 
composed of multiple sub-minima (=tf= 25 and 27 in Tabled] 
- see Sect. I3.1[) . 

4.3. Grand maxima 

4.3.1. Waiting time distribution 

The distribution of the waiting time intervals between sub- 
sequent maxima, listed in Table [H is shown in Fig. [51 with 
the parameters of best-fit approximations shown in Tabled 
row B. The differential distribution (left panel) can be well 
fitted by the power law model. An exponential model gives 
a formally insignificant fit to the distribution. 

The cumulative distribution is shown in the right panel 
and is also close to a power law (see TablelU rowB). The ex- 
ponential model fits short-to-long intervals even better, but 
cannot reproduce the far tail, with three intervals exceed- 
ing 1000 years. Indices for the differential and cumulative 
models are barely consistent with each other (T« r and 
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Table 4. Fitting of power law and exponential models to distributions of the grand minima and maxima occurrence: For 
the differential distribution the value of x 2 is shown together with the corresponding confidence level (in parentheses) 
for 4 degrees of freedom. 





Differential distribution 


Cumulative distribution 


Series 


Power law 


Exponential 


Power law 


Exponential 


A) SN-L 


7 = 1.61 ±0.1 


t = 330 ± 50 yr 


r = 0.95 ± 0.02 


T = 435 ± 15 yr 


min WTD 


X 2 = 0.35 (0.99) 


X 2 = 2.2 (0.7) 






B) SN-L 


7 = 1.36 ±0.1 


T = 430 ± 30 yr 


r = 0.77 ±0.05 


T — 355 ± 20 yr 


max WTD 


X 2 = 0.26 (0.992) 


X 2 = 1.8 (0.79) 






C) SN-S 


7 = 1.82 ±0.06 


r = 250 ± 40 yr 


r = 0.95 ± 0.04 


T — 290 ± 25 yr 


max WTD 


X 2 = 0.22 (0.994) 


X 2 = 6.5 (0.16) 






D) SN-L 


7 = 1.25 ±0.18 


r = 51 ± 5 yr 


r = 1.22 ±0.12 


T = 55 ± 2 yr 


max duration 


X 2 = 1.06 (0.9) 


X 2 = 0.58 (0.97) 






E) SN-S 


7 = 1.5 ±0.6 


t — 50 ± 10 yr 


r = 1.44 ±0.14 


T = 59 ± 6 yr 


max duration 


X 2 = 7.0 (0.14) 


X 2 = 3.3 (0.51) 







i 
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Fig. 7. Histogram of the duration of grand minima. 



y=2.2*x" 

y=2.3*10 3 exp(-x/430) 




- y=41 *x' 
■ y=1 .33*exp(-x/355) 




1000 

Interval length (Years) 



1000 

Interval length (Years) 



Fig. 8. Differential (left panel) and cumulative (right 
panel) distributions of the waiting time between grand max- 
ima according to the SN-L series. 



7 w r+ 1). Accordingly, for the SN-L series we cannot give 
a clear preference to either model. 

The statistics of WTD for the grand maxima using the 
SN-S series is shown in Fig. [9l with the best-fit parame- 
ters listed in Table [4j row C. The power law model sat- 
isfactory fits the differential WTD, while the exponential 
law displays only a poor correspondence to it. The cumula- 
tive WTD (right panel) is nicely fitted by a power law but 
poorly by an exponential. Both models pass the additional 
test (7 w T + 1 vs. T w r). 




~ n . -0.95 

- y=68*x 
y=0.95*exp(-x/290) 



Interval length (Years) 



Interval length (Years) 



Fig. 9. Differential (left panel) and cumulative (right 
panel) distributions of the waiting time between grand max- 
ima according to the SN-S series. 



Therefore we conclude that, although the exponential 
model cannot be totally excluded, the power law model is 
more preferable in describing the WTD of grand maxima. 

4.3.2. Duration of grand maxima 

The distribution of the lengths of maxima in the SN-L se- 
ries is shown in Fig. 1101 with best-fit parameters listed in 
Table [4l row D. The differential distribution (left panel) 
is reasonably fitted by an exponential law but is poorly 
described by a power law. Although both models are seem- 
ingly good in fitting the observed cumulative WTD (right 
panel of Fig. fTT])) . the additional test excludes the power 
law model, since T« r but 7 7^ V ± 1. Therefore, we con- 
clude that the distribution of the lengths of grand maximum 
episodes is close to exponential, as noticed by Solanki et al. 
(2004). 

The differential distribution of the duration of maxima 
for the SN-S series is fitted by none of the two models (see 
Tabled row E). The best-fit parameters for the cumulative 
distribution also favor the exponential distribution (r « T) 
over the power law (T ^ 7 + 1). 

4.4. Quasi-periodicities 

We have also studied possible quasi-periodicities in the rate 
of grand minima/maxima occurrence. We have found that 
the occurrence of grand minima depicts a weak (marginally 
significant) quasi-periodicity of 2000-2400 years, which is a 
well-known period in 14 C da ta (e.g., IDamon fc Sonett 199T1 
|Vasiliev fc Dergachev 2 002). No other periodicities are ob- 
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Duration (Years) 

Duration (Years) 



Fig. 10. Histogram of the duration of grand maxima. 

served in the occurrence rate of grand minima. We have 
found no periodic feature in the occurrence of grand max- 
ima in the SN-L scries, while a marginal hint for a peri- 
odicity of about 1200 years and its harmonics (about 600 
and 400 years - cf. lUsoskin et aT~2"0 04) is found in SN-S 
data. This indicates that the 2400-year periodicity is re- 
lated likely to the clustering of grand minima rather than 
to a long-term "modulation" of solar activity. Therefore, we 
conclude that the occurrence of grand minima and maxima 
is not a result of long-term cyclic variability but is defined 
by stochastic/chaotic processes as discussed in Sect. [5] 

5. Summary of the results 

We have studied the statistics of occurrence of grand min- 
ima and maxima over the last 7-11 millennia. The main 
results can be summarized as follows. 

1. We have presented lists of grand minima (Table QJ and 
maxima (Tables [5] and [3]) , using updated physics-based 
reconstruction of solar activity from 14 C data measured 
by the INTCAL collaboration (jStuiver et al. 1998|) . The 
identification of grand minima is found to be more ro- 
bust to the exact correction of the geomagnetic field 
than the grand maxima. 

2. The occurrence of grand minima and maxima does not 
depict a dominant periodic behaviour. Only a weak ten- 
dency exists for grand minima to cluster with a quasi- 
period of about 2400 years, and no clear periodicities 
are observed in the occurrence of grand maxima. 

3. The waiting time between grand minima depicts a dis- 
tribution which deviates significantly from the exponen- 
tial distribution 0, although the latter cannot be com- 
pletely ruled out because of poor statistics. 

4. The distribution of the duration of grand minima is bi- 
modal, with a dominance of short (30-90 yr) Maunder- 
like minima and a smaller number of long (longer than 
110 yr) Sporer-like minima. 

5. The distribution of the waiting time between grand 
maxima also deviates from an exponential distribution 
(especially for the last 7000 years in the SN-S data se- 
ries), but the latter cannot be completely ruled out. 

6. Lengths of grand maxima correspond to an exponential 
distribution. 

We have tested that the obtained results are robust with 
respect to the uncertainties of the reconstruction. The re- 

3 We confronted the data with a power law WTD, but can 
hardly distinguish it from other non-exponential distributions. 
The most important point here is the deviation from a purely 
exponential WTD. 



suits remain qualitatively the same when varying the pa- 
rameters of Eq. (|A.17[) or when using the sunspot data 
obtained earlier (Solanki et al. 2004, Usoski rTet al. 2006a[) 
based on the old open flux model. 

6. Discussion and Conclusions 

Using the above results we can formulate additional con- 
straints on a dynamo model aiming to describe the long- 
term evolution of solar magnetic activity. 

1. The Sun spends around 3/4 of the time at moderate 
magnetic activity levels (averaged over 10 years). The 
remainder of the time is spent in the state of a grand 
minimum (about 17%) or a grand maximum (9% or 22% 
for the SN-L or SN-S series, respectively). The solar 
activity during modern times corresponds to the grand 
maximum state. 

2. The occurrence of grand minima/maxima is not a result 
of long-term cyclic variations but is defined by stochas- 
tic/chaotic processes. This casts significant doubts on 
attempts of a long-term prediction of solar activity us- 
ing multi-periodic analyses. 

3. The observed waiting time distribution of the occur- 
rence of both grand minima and grand maxima displays 
a deviation from an exponential distribution. A relative 
excess of short and long waiting times indicates that 
the occurrence of these events is not a time independent 
"memoryless" Poisson-like process, but tends to either 
cluster events together or produce long event-free pe- 
riods. Similar waiting time distributions are typical for 
many processes with, e.g. self-organized criticality or 
processes related to accumulation and release of energy 
(see Sect.O. 

4. We distinguish between grand minima of two different 
types: short minima of Maunder type and long min- 
ima of Sporer type fcf. JStuiver fc Braziunas 1989[) . This 
suggests that a grand minimum is a special state of the 
dynamo. Once falling into the grand minimum as a re- 
sult of a stochastic/chaotic but non-Poisson process, the 
dynamo is "trapped" in this state and its behaviour is 
driven by deterministic intrinsic features. 

5. The duration of grand maxima follows an exponential 
distribution, in accord with the earlier finding of Solanki 
et al. (2004). This indicates that leaving a grand maxi- 
mum is a random process, in contrast to the grand min- 
imum case. 

In conclusion, we have presented an analysis of the oc- 
currence of grand minima and maxima of solar activity on 
time scales up to 11,000 years. The results put important 
observational constraints upon the long-term behaviour of 
the solar dynamo. In view of the solar paradigm for the 
magnetic activity of cool stars, we expect these results to be 
applicable also to stellar dynamo models. We note, however, 
that the current results depend on the reliability of the re- 
construction of the sunspot numbers, which in turn depends 
on the reliability of the employed geomagnetic field and 
other factors. This mainly affects the definition of grand 
maxima, while the statistics of grand minima occurrence 
remain fairly robust against these uncertainties. 
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Appendix A: Conversion between the solar open 
magnetic flux and sunspot numbers 

Here we invert the updated model relating the 
sunspot number R to the open magnetic flux F Q 
(Krivova, Balmaceda & Solanki 2007 - referred hence- 
forth as KBS07) to reconstruct the decadal suns pot 
numbers from the open flux (cf. lUsoskin et al. 20*04)) as 
follows. From Eq. (3) of KBS07 one can obtain (henceforth 
(...) denotes 10-year averaging): 



(S) 



(Fa. 



(Fo) 



(F ep h) 



/ dF 
\ dt 

where 

(S) = 

Tta Tte 

From Eq. (1) of KBS07 it follows 

dF ac t /,\ F ac i 

— T— = £act(t , 

at t\ 
where 

111 

T\ T act T t a 



(A.l) 



(A.2) 



(A.3) 



(A.4) 



Hereafter we adopt the parameter values from Table 1 (line 
1) and Sect. 2.1 of KBS07 and express magnetic flux in 
units of 10 14 Wb/month and time in months. Eq. (5) of 
KBS07 takes the form 



Sact(t) w 0.128^(4), 



(A.5) 



Since t± ~ 3 months is much shorter than the cycle length, 
one can assume that F act (t) rs t\ ■ s act (t), and after 10-year 
averaging we obtain: 



(Fact) « 0.38(R). 



(A.6) 



Similarly, since T ep h ~ 0.019 months is very small, one ex- 
pects 



Feph(i) ~ T e ph • £eph(t) , 

where 

(*) = H7 ^act.max.i 

■ sin (f') 



(A.7) 



(A.8) 



Here t' runs over the length of ephemeral region cycle, 
which is longer than a sunspot cycle. With the actual Rq 
data since 1700 ( |Hoyt fc Schatten 1998| ) we have tested 
that the amplitude of a solar cycle is related to the 10-yr 
averaged sunspot value as 



R max = (2.2 ± 0.4) (R). 
Therefore, from (| A.5|) 

tact, max. i ax,i 



(0.28 ± 0.05) (R). 



(A.9) 



(A.10) 



Since e ep h displays long cycles (of the mean duration of 
about 18 years) that partially overlap, the mean value 



of (sin 2 (<')) is numerically found to be 0.74±0.02, and 
Eq. (|A.7|) becomes 



(Feph) « (0.46 ± 0.08) (i?). 



(A.11) 



From the above consideration one obtains that (S) 
(0.0028 ± 0.0001) (R) and 



(R) w (8.5 ± 0.3) • (F ) + (357 ± 13) 



/ dF 
' \ dt 



(A.12) 



In order to evaluate the decadal mean derivative we substi- 
tute the derivative by a slope. 



dF(t)\ A(F) 



dt 



At 



(A.13) 



However, from decadal data we cannot evaluate the value 
of A(F) over a calendar decade: 



A(F) i = (F(t i + 10yT))-(F(t i )), 



(A.14) 



where ti denotes the start year of a decade. Instead we have 
to use the value of 

(F) i+1 - (F)i = (F(t l+1 + 5yr)) - (F{U + 5yr)), (A.15) 

which is displaced in time with respect to the exact defini- 
tion of A(F), From the open flux computed by KBS07 we 
have evaluated that 



dF n 



= (0.73 ±0.06) 



(Fo) 



(Fo) 



(A.16) 



\ dt I ' 120months 

Entering Eq. (|A.16|) into Eq. (|A.12j) . one gets 

(R)i w (8.5±0.3)-(F o ) l + (2.2±0.2)-((F o ) l+1 -(F o ) l ).(A.17) 

We have tested the relation between actual 10-year av- 
eraged GSN and SN computed using Eq. (|A.17[) from the 
open flux (KBS07) for the period 1611-2000. The scatter 
plot shown in Fig. IA.1I displays a good correspondence be- 
tween GSN and SN obtained from Eq. (|A.17|) . The cross 
correlation is r — 0.96, and the difference displays a nearly 
Gaussian distribution with an offset of -0.4 and a = 5.4. 

Thus, we conclude that the conversion between F Q 
and SN decadal data can be done including the influ- 
ence of ephemeral regions in a straightforward manner via 
Eq. (|A.17[) with an uncertainty of a few units in sunspot 
numbers. 
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